Superconductivity in the doped bilayer Hubbard model. 
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We study by the Gutzwiller approximation the melting of the valence bond crystal phase of a 
bilayer Hubbard model at sufficiently large inter-layer hopping. We find that a superconducting 
domain, with order parameter d z 2_ r 2, z being the inter-layer direction and r the intra-layer one, 
is stabilized variationally close to the half-filled non-magnetic Mott insulator. Superconductivity 
exists at half-filling just at the border of the Mott transition and extends away from half-filling into 
a whole region till a critical doping, beyond which it gives way to a normal metal phase. This result 
suggests that superconductivity should be unavoidably met by liquefying a valence bond crystal, at 
least when each layer is an infinite coordination lattice and the Gutzwiller approximation becomes 
exact. Remarkably, this same behavior is well established in the other extreme of two-leg Hubbard 
ladders, showing it might be of quite general validity. 

PACS numbers: 74.20.Mn, 71.30.+h, 71.10.Fd 
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I. INTRODUCTION 

Since its original formulation in the early 60th's, the 
Gutzwiller variational approach^ 2 -! 3 - has proved to be one 
of the simplest yet effective tools to deal with correlated 
electron systems. 

The basic idea of the method is to modify variation- 
ally the weights of local electronic configurations with 
respect to an uncorrelated wavefunction |v?o)i f° r which 
Wick's theorem holds, according to the local interaction 
terms. This is accomplished by means of the variational 
wavefunction: 



(1) 



R 



where Vr is an operator acting on the local Hilbert space 
of the unit cell R. Both the uncorrelated wavefunction 
l^o) an d the operators Vr must be determined varia- 
tionally by minimizing the average energy. In general 
the average energy can be calculated only numerically 4 
but in the limit of infinite coordination lattices, where a 
lot of simplifications intervene 5 - that allow for an explicit 
analytical expressions 7 -^. This is rigorously valid only in 
infinite coordination lattices, nevertheless it is commonly 
used also in finite coordination ones, what is refereed to 
as the Gutzwiller approximation because in a single band 
model it happens to coincide with the approximation in- 
troduced by Gutzwiller himself 3 . 

In spite of its simplicity, many important concepts in 
strongly correlated electron systems have originated from 
Gutzwiller variational calculations or, which is equiv- 
alen1j2ii£, from slave-boson mean-field theory^! 1 .. We 
just mention the famous Brinkmann-Rice scenario^ 2 - of 
the Mott transition. Therefore, even though more rig- 
orous approaches have been developed meanwhile, like 
DMFT 13 or LDA+U 14 , there has been a continuous ef- 
fort towards improving the original Gutzwiller wavefunc- 
tion in finite dimensions 15 , and extending the Gutzwiller 



approximation to account for the exchange interaction 
in multi-orbital models ^ 16 ' 17 ! 18 , for the electron-phonon 
coupling 19 , for interfaces effects 2 ^, and also for more ab- 
initio ingredients 2 ^. A reason for this perseverance is 
that the Gutzwiller wavefunction and approximation are 
so simple and flexible to be adapted to many different 
situations and provide without big numerical efforts rea- 
sonable results. 

In its simplest formulation, Eq. (JTJ) , the form of inter- 
site correlations within the Gutzwiller wavefunction are 
controlled solely by the uncorrelated wavefunction l^o)- 
This aspect should not be problematic if the main interest 
is in the gross features near a Mott transition or when 
a Hartree-Fock Slater determinant gives already a rea- 
sonable description of the actual ground state, which can 
only be improved by applying the operator V . However, 
there are interesting cases where new types of correlations 
may arise near a Mott transition that are not explicitly 
present in the Hamiltonian. A known example are the 
d- wave superconducting fluctuations that are believed to 
emerge in the single band Hubbard model on a square lat- 
tice close to the half-filled antiferromagnetic Mott insula- 
tor—, and which are often invoked to explain high T c su- 
perconductivity. A simple way to justify the emergence of 
superconducting fluctuations is to take the large U limit 
of the Hubbard model, which is known to correspond in 
the low energy sector to the t-J model. Here, the antifer- 
romagnetic exchange J provides an explicit attraction in 
the inter-site singlet channel. This reflects the tendency 
of neighboring sites to form spin singlets, which turns into 
a true antiferromagnetic long-range order at half-filling 
but may mediate superconductivity upon doping. In- 
deed, the Gutzwiller approximation and equivalently the 
slave-boson mean-field theory applied to the t-J model do 
stabilize a <i-wave superconducting phase away from half- 
filling 2 -, a result supported by direct numerical optimiza- 
tions of |*g), with Vr projecting out doubly occupied 
sites and l^o) a d-wave BCS-like wavefunctio n 24 i 25 i 26 i 27 . 
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However, the simplest Gutzwiller approximation in the 
pure Hubbard model away from half-filling does not sta- 
bilize any superconducting phase, just because the on-site 
repulsion U does not couple directly to the d-wave super- 
conducting parameter. A way to improve the wavefunc- 
tion allowing for inter-site spin-singlet correlations could 
be using an enlarged non-primitive unit cell, with "Pn. 
in |T]) acting on a cluster of sites. With this choice the 
wavefunction \^g) breaks explicitly lattice translational 
symmetry, so that one should properly modify the vari- 
ational scheme not to get spurious results, just like any 
other cluster technique^^^L 3 ^ 3 ^. 

Alternatively, one might consider different models that 
are manageable with the simple wavefunction ([I]) and 
which are expected to have a physical behavior similar 
to the one looked for in the Hubbard model on a square 
lattice. One case by now well known is that of two cou- 
pled Hubbard or t-J chains. At half-filling, both models 
are non-magnetic Mott insulator s 34 ' 35 ! 36 . The insulating 
phase is a kind of short-range resonating valence bond 
(RVB) spin-liquid 37 , i.e. a spin-gaped state without any 
symmetry breaking. Actually this state is adiabatically 
connected to the trivial insulator for very large inter- 
chain coupling, which is a collection of inter-chain dimers, 
what can be denoted as a valence bond (VB) crystal. 
Away from half-filling, dominant superconducting fluc- 
tuations aris o 36 ' 38 ' 39 , with a two-chain analog of a two- 
dimensional d-wave symmetry. The emergence of strong 
superconducting fluctuations appears here as the natural 
fate of doping the half-filled VB Mott insulator—, real- 
izing in one-dimension the RVB superconductivity sce- 
nario proposed by Anderson 4 ^ in the early days after the 
discovery of high T c superconductivity. An immediate 
question that arises is whether the above one-dimensional 
behavior survives in higher-dimensions, namely how ro- 
bust is the two-chain RVB scenario upon increasing di- 
mensionality. This is actually the content of the present 
work. 

As a matter of fact, this question has been already 
addressed several times in connection with high T c su- 
perconductors, specifically analyzing a bilayer Hubbard 
model by various techniques, including quantum Monte 
Carlci 1 -! 4 ^ 3 .! 4 ^ (QMC) and DMFTi 6 -^. In section HI 
we shall discuss more in detail these early works while 
introducing the model. More recently, the same problem 
has been studied at half-filling by an improved Gutzwiller 
approximation 1 ^, which we present in section lHll together 
with a further improvement that we use here to extend 
that analysis away from half-filling. The results are pre- 
sented in section [TV] while section IVl is devoted to con- 
cluding remarks. 



II. THE MODEL 

Throughout this work we shall be interested in a bi- 
layer Hubbard model described by the following Hamil- 



tonian 

2 

ft = - C R,i<r C R',i<7 + H - c - 

RR' i=l a 

+|£I>*m- 1 ) 2 

R i = l 

J2 ( C R,la C R,2 CT + #-c) 

R a 

2 

= e(k ) c L<t c m<t 

kcr i— 1 

+yE -i) a 

R i=l 
Via 

= 7~tho P + 7~Lu +7~L±, (2) 

where t± > 0, cj^ io . and c R ia create and annihilate, re- 
spectively, an electron at site R in plane i = 1, 2 with spin 

0") n R,i — So- c r io- c R ia 1S * ne l° ca l occupation on layer 
i, and U is the Hubbard repulsion on each lattice site. In 
order to study the doped system it is more convenient to 
work in the grand-canonical ensemble adding a chemical 
potential term — ^Sr^R^ to the model Hamiltonian 
([2]). The particle number is then controlled by tuning fi. 
In Eq. ^) ia creates an electron in layer i and spin a 
with momentum k, and e(k) E [-D, D] is the intra-layer 
dispersion in momentum space, where D is half the band- 
width that will be our unit of energy. The non-interacting 
part of the Hamiltonian is better rewritten introducing 
the bonding (e) and antibonding (o) combinations 

c f = — ( c f + c f 

L k,ecr ~ ^ V k,l<T ^ c k,2a J ' 
c f - — ( c f - c f ) 
through which 

n hop + H± = J2 £ e «( k ) c U c k,«a. ( 3 ) 

Via a—e^o 

where e e (k) = e(k) — t± € [— D — t±, D — t±] and 
e Q (k) = e(k) + t± e [— D + t±,D + t±] are, respectively, 
the bonding and antibonding band dispersions. 

If U — and the density is one electron per site, half- 
filling, the model describes a metal until the two bands 
overlap, i.e. t± < D, and a band insulator otherwise. 

For U ^> D + t±, the model becomes equivalent to two 
Heisenberg planes coupled to each other by an inter-plane 
antifcrromagnetic exchange J± = At\/U. If each plane 
is a square lattice with only nearest neighbor hopping 
t, hence D = 4i, each Heisenberg model is character- 
ized by a nearest neighbor antiferromagnetic exchange 
J = 4t 2 /U. This model has been studied in detail by 
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quantum Monte Carlo 4 and it is known to have a 
quantum critical point that separates a a Neel antifer- 
romagnet, for J± < 2.5520 J, from a gaped spin-liquid 
phase, for larger J±. The latter can be interpreted as 
a kind of VB crystal, each bond being an inter-layer 
singlet, adiabatically connected to the band insulator at 
U = 0. In terms of the hopping parameters of the origi- 
nal Hubbard bilayer, the critical point should correspond 
to {t±/t) c = v / 2.5220 ~ 1.5881. This value is in good 
agreement with direct QMC simulations of the Hubbard 
bilaye r 42 ' 45 , which find (t±/t) c ~ 1.5 to 2. According to 
these results, when 1.6 < (t±/t) < 4 one could start at 
U = with a metallic phase, and, upon increasing U, find 
a direct transition into the VB Mott insulator. However, 
the story must become more complicated if the U = 
Fermi surface at half-filling has nesting at the edge of 
the Brillouin zone, as it happens for a square lattice with 
only nearest neighbor hopping. In this case, the U = 
and t± < At = D metal has a Stoner instability towards 
Neel antiferromagnetism for arbitrary small U, so that 
it is a priori not obvious that one could find any direct 
metal to VB Mott insulator transition. In reality, both 
cluster DMFT 47 and QMC simulations find evidence that 
such a transition does exist. Nevertheless, one may al- 
ways bypass this problem assuming that the intra-layer 
hopping is such as not to lead to any nesting, the latter 
being more an accident than the rule in realistic systems. 
In this case, which we will implicitly assume hereafter, it 
is safe to believe that a direct transition at half-filling 
from a metal to a VB Mott insulator does exist. 

Within this scenario, the melting of the VB crystal 
into a metallic phase can therefore occur either by dop- 
ing away from half-filling but also upon decreasing U 
below the Mott transition, still keeping half-filled den- 
sity. In the latter recent studyi 8 - has shown that, 
within the Gutzwiller approximation, the VB crystal first 
turns into a superconducting phase that eventually gives 
way to a normal metal upon further decreasing U . This 
finding supports the RVB superconductivity scenario^ - 
and shows that the one-dimensional behavior persists in 
higher dimensions. It also agrees with the indication 
of an enhanced pairing susceptibility obtained in earlier 
studies by QM C 43 ' 48 . However the lowest temperatures 
attainable so far by QMC are still above the eventual 
superconducting critical temperature, so that the exis- 
tence of a true superconducting phase at half-filling is 
numerically still an open issue. DMFT calculations, that 
could in principle be carried out at zero temperature, 
was performe d 46 ' 47 but did not search explicitly for any 
superconducting phase. 

Away from half-filling, QMC indications of enhanced 
pairing fluctuations are more convincin g 43 ' 45 , although 
the existence of a superconducting phase at low temper- 
ature is still uncertain^ 3 -. This makes it worth addressing 
this issue by the Gutzwiller approximation, which is not 
as rigorous as QMC but at least can provide results at 
zero temperature. 



III. THE METHOD 

In order to study the bilayer Hubbard model (JSj) away 
from half-filling we adopt the Gutzwiller approximation 
scheme developed in Refs.ll8lto deal with the same model 
at half-filling. The variational wavefunction that we use 
has the form as in Eq. (TTJ) where 

1. Vr acts on the full Hilbert space that includes site 
R in layer 1 and site R in layer 2; 

2. I^o) is allowed to be a BCS-wavefunction with sin- 
glet order parameter in the channel c R ^4 2 j, + 

J J " 

L R,2T L R,11- 

The most general expression for Pr is: 

V R = M R )rir 2 l r i,R)(r 2 ,R|, (4) 
ri,r 2 

where each state |rj,R) denotes a local two-site elec- 
tronic configuration, and the matrix A(R) has to be vari- 
ationally determined. Average values of operators on 
the wavefunction (fTJ can be analytically computed in 
infinite-coordination lattices provided the following con- 
straints are satisfied by Prt ' 17 ' 18 : 

(*o|P r 7>r|*o) = 1, (5) 

(*o|7> r 7> r Cr|* ) = <*o|Cr|* ), (6) 

where Cr is the local single-particle density-matrix op- 
erator with elements 4a c R/3 and 4a4fl> a and (3 
labeling single-particle states (both layer and spin in- 
dices) and c Ra (c RQ ) creating(annihilating) an electron 
in state a at site R. Expectation values of local operators 
are then computed straightforwardly-^, as (V^OrV )o = 
(7 5 r O r 7 : 'r)o (hereafter (...)o denotes averages on the 
uncorrelated wavefunction \^q), which can be easily com- 
puted by means of Wick's theorem). When calculating 
the average of the inter-site density matrix, one finds that 
the physical single- fermion operator acting on \^>g) is ef- 
fectively replaced by a renormalized one acting on |\&o) 
according to: 

4,a -> E 4/3 + E QWafi C R,/3' ( ? ) 

P 

where the renormalization matrices R and Q are deter- 
mined by inverting the following set of equations 18 : 

(^R c R,a^R c R,/3>0 = E fl ( R )«7 ( C R )7 C R,/3 >0 

7 

+ E <2( R )«7 ( C R,7 C R,/3 >0> (8) 

7 

<^r4,^r4 i/3 >o = E fi ( R )«7 <4, 7 4,/3 >o 

7 

+ E <2( R )«7( C R, 7 4,/3>0- (9) 

7 
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In Ref. [l8|, a different notation was used for the ma- 
trices R and Q, namely R — y/Z and Q = \/~K. In 
order not to generate any confusion with the definition 
of a square root of a matrix, and also for keeping more 
explicit the connection with slave-boson mean field the- 
oryii, we have preferred here to use R and Q. Despite 
the considerable simplification introduced by the infinite- 
coordination limit, the variational problem remains still 
a difficult task to deal with because of the large size of 
the local Hilbert space, which contains 16 states so that 
A spans in principle 16 x 16 matrices. 

A further simplification can be achieved with a proper 
choice of the basis set spanning the local Hilbert space. 
This can be done, for instance, by using from the begin- 
ning the natural basis, i.e., the single-particle basis which 
diagonalizes the variational density matrix (Cr,)ct^. An 
alternative and more efficient approach consists^ in 
defining the local operator Vr in a mixed-basis repre- 
sentation, namely expressing |ri,R) = |T, R) in Eq. ((3]) 
in the original basis defined by the model Hamiltonian 
and assuming that (r2,R| = ({n Q },R| are Fock states 
in the natural basis, identified by the occupation numbers 
h a = 0, 1. With this choice, one can use as variational 
parameters just the eigenvalues of the density matrix, be- 
cause the unitary transformation that relates the natural- 
basis operators d R to the original ones c R needs not 
to be known explicitly. This simplifies considerably all 
calculations. In the mixed original- natural basis repre- 
sentation one introduces a new matrix 



operators and the renormalization factors acquire very 
simple expressions: 



(R) = A(R) /P°(R), 



(10) 



where A(R) is the variational matrix in the mixed- 
basis representation and P°(R) is the uncorrelated 
occupation-probability matrix, with elements 

P^ a}M (K) ee ( \{m a }, R)({n Q },R| ) 

being 

P°(R) { ^ } = [J (n°(R) Q ) fiQ (1 - r^R)*) 1 ^" , (11) 

a 

and n°(R) Q the eigenvalues of the density matrix (to be 
variationally determined). In terms of <j>, the constraints 
to be imposed on the Gutzwiller wavefunction can be 
recast asM: 

Tr(^) = 1, (12) 

Tr (ftftdidp) = 6 a , n° a , (13) 

Tr(>V44) = °> ( 14 ) 

where, to simplify notations, we dropped the site-label 
R, and we also introduced matrix representations of the 
fermionic single-particle operators. The average of local 



[V^OV) Q 

Raf3 
Qa/3 



Tr (tfoA , (15) 
1 TrjVct^ Y (16) 

1 TrLuUdl). (17) 



Through Eqs. (fl6|) and (fT7|) single-particle fermionic op- 
erators are automatically mapped into renormalized op- 
erators in the natural basis, and Eq. ([7]) is replaced by: 



(18) 



^ Ra/3 dp + Qa/3 dp. 



Note the presence of the latter term in the rhs of Eq. (| 18(1 . 
which makes it possible that a creation operator in the 
original representation turns into an annihilation oper- 
ator in the natural one. Its existence is a direct conse- 
quence of allowing |\& ) to span also BCS-like wavefunc- 
tions and/or Vr to couple states with different particle 
numbers. Should \^f Q ) describe a normal metal and Vr 
be diagonal in the particle number, Q a p would be strictly 
zero, as was the case in Ref.l49l Therefore Eqs. (fTrJl) . (TIT)) 
and HU) extend Eqs. (A10) and (A18) of Ref. El to the 
more general case in which superconductivity is allowed. 

Practically, it is convenient^ to generate variational 
matrices <j> that directly satisfy Eqs. (fl"2|) - (fl"4"|) hence uni- 
vocally determine the parameters rfi a , and only after im- 
pose, by proper Lagrange multipliers, that the uncor- 
related l^o) has an average local density matrix with 
eigenvalues n a . We end mentioning that the elements 

0r,{B} = Ar{n} y/Pj^ of the matrix (|10[) correspond to 
the slave-boson saddle-point values within the mean-field 
scheme recently introduced by Lechermann and cowork- 
ers^. 

It may happen that, in spite of all the above simplifica- 
tions, the variational space thus generated is still unnec- 
essarily large. For instance, if one looks for a variational 
wavefunction that preserves particle number, all the el- 
ements of A connecting subspaces of the local Hilbert 
space with different particle numbers should be identi- 
cally zero. Therefore it would be desirable to specialize 
the general procedure sketched above in such a way that 
symmetries can be built in the variational wavefunction 
from the onset. In general, given a symmetry group Go 
that one would like to enforce, we must require, in addi- 
tion to (O and (HI), that 



[Pr, Go] = 0, 



(19) 



However, in the mixed representation there may be some 
symmetry operations that can not be defined without an 
explicit knowledge of the natural basis in terms of the 
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original one, which would make the whole method much 
less efficient. If one decides not to implement these sym- 
metries, but only those, symmetry group G C Go, that 
commute with the most general unitary transformation 
U connecting original and natural basis, i.e. 

[U,G] = 0, 

compatibly with the variational ansatz, the above de- 
scribed variational method can be still used with the fol- 
lowing modification. 

Let us assume this case and define a unitary operator 
V that transforms the Fock states in the original basis 
into states that decompose the local Hilbert space in ir- 
reducible representations of the group G. We define G 
the representation of G in such a basis. Because of our 
choice of the subgroup G, V does right the same job even 
in the natural basis, although this is unknown. Since 
the trace is invariant under unitary transformations, all 
formulas from (|12|1 to (fl~7|) remain the same even if the 
variational matrix cj> and the matrix representation of the 
single fermion operators are defined in the states of the 
irreducible representation, both in the original and nat- 
ural basis, with the additional symmetry constraint 



,G] =0, 



(20) 



which follows from ([19]) . We note that the matrix rep- 
resentation of a single-fermion operator in these states 
is readily obtained once V is known, and is trivially the 
same for both original and natural operators. Therefore 
it is sufficient to create and store it at the beginning of 
any calculation. 

As an example, which is directly pertinent to this work, 
let us consider G the group of spin SU(2) transforma- 
tions. In this case an irreducible representation is read- 
ily obtained and consists of states with fixed total spin 
S and its z-projection S z , of the general form |T, S, S z ), 
where T serves as an additional label to distinguish be- 
tween different states with same S and S z in the original 
representation. The operator V is thus the unitary trans- 
formation that connects Fock states in the original basis, 
|{n a }) to the states |r, S, S z ), 



controlled by a chemical potential term -C^Ri^R-,' 
that we add to the Hamiltonian @. We search for a 
variational solution that allows for singlet superconduc- 
tivity and doesn't break spin-S'?7(2) symmetry. In this 
case, the unitary transformation that would connect orig- 
inal to natural bases, would leave spin SU (2) generators 
invariant, so that the method described in the previous 
section is applicable. 

We already said that the Gutzwiller operator Vr in 
(HJ acts on the whole local Hilbert space of two sites, R 
in layer 1 and R in layer 2. The variational energy to 
be minimized is then the sum of two terms. One is the 
contribution of the local, same R but both layers, terms, 
which, according to the results of the previous section, 
reads: 



Eioc = Tr 



6+ \Hu + H ± 



Ri 



(22) 



where all operators are meant to be matrices in the lo- 
cal representation invariant under SU (2) symmetry. The 
other contribution to the total energy is the intra-layer 
hopping Ehop- This can be shown to coincide with 
the ground-state energy of a variational single-particle 
Hamiltoniani^: 



n *ho P = X^k^k^k, 



(23) 



where t/>£ = (^k,it'^k,2T>^-kij.)^-k,2j.) 1S ^ ne Nambu 
spinor in momentum space and Tk a 4 X 4 matrix in 
the natural basis that depends explicitly on momentum 
and on some Lagrange multipliers. These are included 
to enforces that the average of the single particle density 
matrix on the ground state - to be identified with l^o) 
in (JTJ) — is diagonal with matrix elements satisfying 



(*oi4 



V : \{n a }) 



\T,S,S Z 



The matrix Tk has the general expression: 

f ( e(k)Z + fj e(k)A + S \ . . 

k 1 e(k)At+# - e (k)i*-j7* J ' { ' 



We use the same V to generate from the Fock states in the 
natural basis, which we remind is and remains unknown, 
the states \T, S, S z ). It follows that, in order to preserve 
full spin SU (2) symmetry, 

^T,S,S Z ;T,S,S Z — $SS ^S Z S Z ^r,S;f,S' (21) 

is a block matrix. 



IV. RESULTS 

Let us turn now to study the bilayer Hubbard model 
@ away from half-filling. As mentioned, the filling is 



where the 2x2 matrices fj and 6 are the aforementioned 
Lagrange multipliers, while Z and A have elements (la- 
belled by j, I = 1,2, the layer indices) 

2 

y-w = E ( R i,j R t,i-Qi,iQij)> ( 25 ) 



2 



i=i 



(26) 



We solved numerically the variational problem as- 
suming for simplicity a flat density of states with half- 
bandwidth D (we do not expect the results to change 
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qualitatively by adopting a more realistic density of 
states). In order to compare with the half- filling results 
reported in Ref. HH, we fixed the value of the intra-dimer 
hopping t±/D = 0.5 and solved the variational problem 
for different values oiU/D and n/D. Note that this value 
in the case of a square lattice with nearest neighbor hop- 
ping t corresponds to ij_ = 2t, above the critical value 
for the stability at large U of the VB Mott insulator—. 

At half- filling, n/D = and we recover all results of 
Ref. Specifically, we find a first order metal to VB 
insulator transition. In the metallic phase just before the 
transition, singlet superconductivity emerges. In Fig. [T] 
we show as function of U/D the behavior of the inter- 
layer, Aj_ and intra-layer, A||, superconducting order pa- 
rameters, defined as 

A± = (*g| 4,^4,2; + 4,2t4, u I*g>' (27) 
A|| = 4,iT4^ + 4', iT 4 >i ;l*G>, (28) 

where R and R' are nearest neighbor sites on layer 
i = 1,2. We find that, near the first order transition 
that we think identifies the actual Mott transition, both 
order parameters are finite and have opposite sign, the 
so-called d z 2_ r 2 symmetry known to be dominant in the 
two-chain model^, and which QMC simulations^^ in- 
dicate as the leading pairing instability. The variational 
energy that we obtain appears to be slightly lower than 
that found in Ref. uM, as one could have expected due to 
the larger number of variational parameters. Nonethe- 
less, the critical U c at the Mott transition is only slightly 
reduced to U c /D ~ 2.02 for tj_/D = 0.5. We note that 
the phase at U > U c , that we believe is Mott insulat- 
ing, still shows a finite superconducting order parameter 
that dies out upon increasing U. As discussed in[l8|, we 
think this might be a spurious result of our variational 
approach that lacks intersite charge correlations crucial 
in stabilizing a genuine Mott insulating phased. 

We study finite hole doping by varying [i/D < at dif- 
ferent values of U/D. Before discussing the variational 
results, we briefly sketch the behavior of the doped non- 
interacting system, U/D = 0. The inter-layer coupling 
gives rise to bonding and antibonding bands, see Eq. ([3]). 
With the chosen value of t± = 0.5D, these bands overlap 
at half-filling and the system displays a metallic behavior. 
When the chemical potential is lowered, holes are injected 
into the system inducing a depletion of both bands un- 
til, at a given value of the chemical potential, the upper 
(antibonding) band empties. For the chosen t± and for 
a flat density of states the complete depletion of the an- 
tibonding band happens at fi — 0.5D, corresponding to 
quarter filling n — 1. As a consequence, both the intra- 
layer (Ehop) and inter-layer (E±) hopping contributions 
display a discontinuity in their first derivatives at quarter 
filling, signaling that the antibonding band is no longer 
contributing. The total energy however remains smooth 
for any value of fi (or equivalently n), as it should. When 
U/D 0, the behavior that we find depends crucially if 
U is smaller or greater than U c , namely if the half-filled 
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FIG. 1: (Color online) Inter- plane (blue circles) and, with re- 
versed sign, intra-plane (red triangles) superconducting order 
parameters at half-filling as function of U/D. The vertical 
line indicates the first order transition that we think identi- 
fies the on-set of Mott insulating behavior. Inset shows the 
variational energy in units of D 
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FIG. 2: (Color online) Average density n summed over both 
layer as a function of the chemical potential /i < for selected 
values of interaction U/D. 

state is a metal or an insulator. 



A. Doping the metal at U < U c 

As long as U < U c , any change of /i induces a continu- 
ous change in the total particle number; a finite com- 
pressibility signal of a metallic behavior, as shown in 
Fig. [2l Alike the uncorrelated case, a cusp appears in 
the evolution of n at quarter-filling, that we explain seem- 
ingly as the depletion of the antibonding band. Indeed, 
when U < U c , the metallic solution evolves just like the 
non-interacting case. The main effect of interaction is to 
slightly reduce inter- and intra-layer hopping contribu- 
tions with respect to their uncorrelated counterparts, as 
shown in Fig. [3] where we plot the different contributions 
Ehop, E± and Ejj to the variational energy. The intra- 
layer hopping contribution Eh op diminishes in absolute 
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FIG. 3: (Color online) Left panel: The different contributions 
to the variational energy as a function of doping for U/ D = 1 
and per lattice site R, i.e. summed over both layer. As a ref- 
erence, the behavior of non-interacting inter- and intra-layer 
hopping contributions is plotted (dotted lines). In the inset 
the total energy E var (n) = E var + fin is shown: despite the 
cusp observed in the hopping contributions, the evolution of 
E var (n) is smooth. Right panel: Occupation of the varia- 
tional lower and upper bands as function of n. Dotted lines 
represent average occupation of even and odd orbitals. 



value with increasing doping because of the depletion of 
the bands, as it occurs in the non-interacting system; 
at quarter-filling it displays a cusp and correspondingly 
the inter-layer hopping E± starts to rapidly decrease, the 
effects of U being more and more negligible as the low- 
density regime is approached. In the right panel of Fig. [3] 
we show the occupancies n® and of the variational 
lower and upper bands, respectively, which are obtained 
by diagonalizing the associated variational Hamiltonian, 
Eq. (|23p . and actually coincide with the eigenvalues of 
the single-particle density matrix. As in the uncorre- 
cted system, the occupancy of the upper band vanishes 
at quarter filling. We stress the fact that in the present 
approach these states are variationally determined and 
may not correspond to the even and odd combinations 
of the original operators. However, as long as U < U c , 
we find that the average values of bonding and antibond- 
ing band occupancies, n e and n Q , almost coincide with, 
respectively, nf and n°. 

Concerning superconductivity, we find that the inter- 
layer order parameter, Eq. (|27p . is extremely small, prac- 
tically zero within our numerical precision, see Fig. 3] 
The intra-layer order parameter strictly follows the inter- 
layer one, hence is also zero. 



B. Doping the VB Mott insulator at U > U c 

When U > U c , i.e. when the half-filled system is insu- 
lating, the particle number remains stuck to its half-filled 
value n = 2 until \fi\ < \fi*\ m (U - U c )/2. This simply 
follows from the existence of the Mott gap at half-filling. 
Upon doping, i.e. when \fi\ > a metallic behavior is 
clearly found. However, within our numerical precision 
we can not establish whether the evolution from the insu- 
lator to the metal occurs smoothly (yet with a diverging 



compressibility) or through a weak first-order transition. 
Till the largest value of U we considered, we could not 
find any appreciable discontinuity in the evolution of n 
at large doping, unlike for U < U c where a cusp is ob- 
served at quarter filling. In addition, contrary to the 
case U < U c , here we find a clear superconducting signal 
between half and quarter filling, see e.g. the behavior of 
Aj_, Eq.[23 shown in Fig.0J We note that Aj_ has a non- 
monotonous behavior, first increases quite rapidly with 
U and for larger values decreases. Like at half-filling, a 
finite Aj_ produces through Eq. (fT5|) also a finite intra- 
layer A||, Eq.f251 not shown here, which happens to have 
opposite sign. 



Let us now consider in detail the energetic balance for 
U > U c and its differences with respect to U < U c . At 
very large U (not shown), as holes are injected into the 
system, both intra- and inter-layer hopping contributions 
first increase in absolute value, then saturate around ap- 
proximatively quarter-filling, and eventually decrease as 
the low-density regime is attained, as expected when ap- 
proaching the bottom of the variational bands. In other 
words, the behavior at large U between half- and quarter- 
filling is quite different from the non-interacting case, 
while becomes quite similar below. This points to a 
very different influence of a strong interaction close to 
half-filling and far away from it and, indirectly, empha- 
sizes the role of the superconductivity that we find for 
2 > n > 1. For U > U c , i.e. closer to the half-filled 
metal-insulator transition, the picture is slightly differ- 
ent, as shown in Fig. O for U/D = 3. To begin with, 
at small dopings the system gains in intra-layer hopping 
energy while the inter-layer one seems to be slightly re- 
duced. Remarkably, even if the total energy is, within our 
numerical accuracy, a smooth function of n, both hop- 
ping contributions display a discontinuity at n/D ~ 1.28, 
which corresponds to a local density of n « 1.27. Here 
the occupation of the upper variational band goes to zero 
(cfr. right panel of Fig. [5]), even though nothing simi- 
lar occurs in the occupation of the physical antibonding 
band. At this filling fraction, the inter-layer hopping en- 
ergy gain has an upward jump, contrary to the intra-layer 
one, even though further doping leads to a reduction of 
both. A drop in the amplitude of the superconducting 
order parameter Aj_ is also found at this point. Further 
doping diminishes Aj_, which vanishes approximatively 
at quarter filling. A similar feature is observed in an- 
other quantity. Indeed, just like nf and n° may not 
correspond to the occupation of the bonding and anti- 



bonding bands, n° = rtj 1 



which is the average den- 



sity of the BCS-likc variational wavefunction, may differ 
from the physical one. In the inset of Fig. [4] we show 
their difference for U/D = 3. We observe that they ac- 
tually deviate when superconductivity is found and their 
difference jumps down abruptly for n < 1.27. 
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FIG. 4: (Color online) Superconducting inter-layer order pa- 
rameter Ax for different U / Ds. In the inset we plot the differ- 
ence between the local densities of the BCS variational wave- 
function |\&o) and of the actual one I^g), at U/D — 3. 
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FIG. 5: (Color online) Left panel: Contributions to the varia- 
tional energy as function of doping for U/D = 3. In the inset 
the total energy E var {n) = E var + fin is shown: despite the 
discontinuities observed in the hopping contributions, the evo- 
lution of E(n) is, to our numerical accuracy, smooth. Right 
panel: Occupancies of lower, and upper, n°, variational 
bands as function of n. Dotted lines represent average oc- 
cupation of even, n e , and odd, n , orbitals. Note that the 
insulating phase at half-filling is identified by the lower band 
fully occupied and the upper one empty. The latter empties 
again for doping 2 — n > 0.73. 



bilayer Hubbard model. We have considered a value of 
the inter-layer hopping t± such that, at half-filling, the 
model should undergo a direct transition at U — U c from 
a metal to a non-magnetic Mott insulator, a valence bond 
crystal consisting of inter-layer dimers. This choice offers 
the opportunity to study how a valence bond crystal liq- 
uefies either by reducing the Coulomb repulsion keeping 
the density fixed at one electron per site, or by adding 
mobile holes. The melting upon decreasing U was already 
shown 1 — to lead to a superconducting phase intruding be- 
tween the valence bond insulator at large U > U c and the 
normal metal at weak U <^iU c . Here we show that super- 
conductivity arises also upon melting the valence bond 
crystal by doping. In other words, the superconducting 
dome that exists at half-filling close to U c extends into a 
whole region at finite doping. The maximum supercon- 
ducting signal is found at 20% doping, and beyond that 
it smoothly diminishes, disappearing roughly at quarter 
filling within our choice of parameters. These results are 
appealing as they show that the well established behav- 
ior of a two-leg Hubbard ladde r 18 i 36 ' 37 i 39 seems to survive 
in higher dimensions, actually in the infinite-dimension 
limit where our Gutzwiller approximation becomes ex- 
act. It is obvious that, in spite of all improvements of the 
Gutzwiller variational approach, to which we contribute a 
bit with this work, this method remains variational hence 
not exact. Therefore it is still under question if supercon- 
ductivity indeed arises by metallizing the valence bond 
Mott insulating phase of a Hubbard bilayer, which we 
believe is an important issue of broader interest than the 
simple bilayer model we have investigated^. There are 
actually quantum Monte Carlo simulation :/ 2 ' 43 ' 45 ' 48 that 
partially support our results as they show a pronounced 
enhancement of superconducting fluctuations close to the 
half-filled Mott insulator. However a true superconduct- 
ing phase is still unaccessible to the lowest temperatures 
that can be reached by quantum Monte Carlo. On the 
other hand, dynamical mean field calculations, that can 
access zero temperature phases, did not so far looked for 
superconductivit y 46 ' 47 . Therefore we think it would be 
worth pursuing further this issue. 



V. CONCLUSIONS 



In this work we have studied by means of an extension 
of the Gutzwiller approximation the effect of doping a 
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